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Preface 


The  purpose  of  this  study  project  was  to  develop  a  simple  method 
of  predicting  fallout  from  a  low  yield,  nuclear  airburst  that  would  be 
useful  to  a  tactical  ground  commander.  Simple  calculation  methods 
exist  for  tactical  use;  however,  none  are  utilized  for  airbursts.  A 
computer  model  was  developed  based  on  low  yield  test  data.  The  method 
produces  maximum  dose  rate  and  infinite  dose  curves  enabling  a  commander 
to  quickly  identify  radiation  hazard  areas. 

I  thank  LCDR  James  H.  Gogolin  and  Capt  Arthur  T.  Hopkins  for  their 
contributions  to  the  lively  art  of  fallout  conversation  and  providing 
constructive  criticism  during  this  study.  I  particularly  want  to  express 
my  appreciation  and  thanks  to  Dr.  Charles  J,  Bridgman  for  his  guidance 
and  sage  advice  in  insuring  that  the  cart  always  stayed  behind  the  horse. 
Lastly,  I  want  to  thank  my  wife,  Elsa,  who  always  made  sure  that  there 
was  a  light  still  flickering  at  the  end  of  the  tunnel. 
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Abstract 
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A  method  was  developed  that  enables  a  tactical  ground  commander  to 
predict  gamma  radiation  dose  rates  and  infinite  doses  produced  by  the 
rain  scavenging  of  low-yield  nuclear  airburst  clouds.  A  ground  activity 
distribution  per  unit  area  at  time  t,  A(x,y,t)  ,  is  computed  using  a 
distribution  of  activity  in  the  cloud  per  meter  of  altitude,  A(z,t) 

To  find  the  maximum  activity  grounded,  it  is  assumed  100  percent  of  the 
cloud  activity  is  instantaneously  deposited  on  the  ground  by  the  mechanism 
of  rain  scavenging.  This  maximum  A(x,y,t)  is  then  converted  to  a  maximum 
dose  rate,  D(x,y,t),  from  which  maximum  infinite  doses  are  computed. 
Maximum  dose  rate  and  infinite  dose  curves  vs  cloud  washout  time  (after 
cloud  stabilization)  are  presented  for  10  weapon  yields  ranging  from 
1.0  kiloton  to  0.1  kiloton.  It  is  shown  that  the  radiation  hazard  levels 
are  insignificant  tactical  threats  at  times  greater  than  36  hours  after 
cloud  stabilization.  .. 


A  RAIN  SCAVENGING  MODEL  FOR  PREDICTING 
LOW  YIELD  AIRBURST  WEAPON  FALLOUT  FOR 
OPERATIONAL  TYPE  STUDIES 

I .  Introduction 

The  ability  to  predict  fallout  from  the  atmospheric  detonation  of 
low  yield  nuclear  weapons  plays  an  important  role  in  tactical  planning 
and  operations.  A  tactical  ground  commander  must  have  the  capabilif  o 
predict  radiological  contamination  hazard  areas  resulting  from  the 
friendly  or  unfriendly  employment  of  nuclear  weapons.  These  hazard  a  o 
are  capable  of  producing  mass  casualties.  Since  actions  must  be  taken 
to  minimize  the  casualty  producing  effects  of  these  hazard  areas  to 
friendly  forces,  the  tactical  ground  commander  must  be  able  to  predict 
the  location  and  intensity  of  the  hazard. 

Currently  a  method  exists  that  allows  ground  commanders  to  predict 
fallout  hazard  areas  for  all  types  of  nuclear  weapon  bursts  except  an 
airburst  (Ref  5).  When  an  airburst  is  employed  no  tactically  signifi¬ 
cant  fallout  will  occur  because  airburst  fallout  particles  are  too  small 
to  be  deposited  locally  on  the  battlefield.  However,  if  a  rain  (or  snow) 
cloud  interacts  with  the  nuclear  fallout  cloud  (scavenging  the  cloud), 
then  tactically  significant  fallout  could  result.  A  simple  method  use¬ 
ful  to  a  tactical  ground  commander  for  predicting  low  yield  weapon  fall¬ 
out  due  to  rain  scavenging  does  not  currently  exist  (Ref  6:B-12). 


This  thesis  develops  a  simple  rain  scavenging  fallout  model  for 
the  tactical  commander.  For  the  purpose  of  this  thesis,  low  yields 
will  be  defined  as  being  less  than  or  equal  to  1  kiloton.  Section  II 
provides  an  introduction  to  the  two  principle  types  of  fallout  models 
currently  in  use.  Section  III  explains  the  concept  of  the  stabilized 
cloud  and  describes  a  low  yield  cloud  model  based  on  empirical  weapons 
test  data.  Section  IV  describes  the  rain  scavenging  model  for  airbursts 
Section  V  contains  the  results  of  this  model  and  presents  a  method  which 
a  tactical  ground  commander  can  use  in  assessing  the  radiological  hazard 


to  his  troops. 


This  section  describes  the  two  types  of  fallout  models  presently 
in  use.  The  first  type  is  the  numerical  model  represented  by  the 
Department  of  Defense  Land  Fallout  Interpretive  Code  (DELFIC)  (Ref  12). 
The  second  type  is  the  analytical  model  represented  by  both  the  Pentagon 
Weapon  Systems  Evaluation  Group  (WSEG)  and  Air  Force  Institute  of  Tech¬ 
nology  (AFIT)  codes  (Refs  14,  3). 

Numerical  Model  (DELFIC) 

DELFIC  is  a  full  physics  computer  code  developed  for  use  as  a 
research  tool  (Ref  12:2).  DELFIC  models  the  fallout  process  by  comput¬ 
ing  the  space-time  history  of  fallout  particles  as  they  are  transported 
through  the  atmosphere  and  down  toward  the  earth's  surface.  These  fall¬ 
out  particles  are  represented  by  panca)te  shaped  wafers  containing  distinct 
particle  size  groups.  The  fallout  pattern  on  the  ground  at  a  given  time 
is  computed  by  superpositioning  the  ground  positions  of  those  wafers 
that  have  reached  the  ground  in  that  given  time.  Because  DELFIC  must 
carry  out  a  large  number  of  calculations,  it  is  a  slow  running  code  and 
very  expensive  to  use.  To  alleviate  this  problem,  faster  running  ana¬ 
lytical  codes  have  evolved. 

Analytical  models  (WSEG  and  AFIT) 

The  WSEG  and  AFIT  codes,  unlilte  DELFIC,  model  the  fallout  particle 
deposition  by  smearing  the  stabilized  radioactive  cloud  on  the  ground 
as  it  falls.  The  resulting  fallout  pattern  on  the  ground  is  illustrated 


by  the  shaded  area  in  Figure  1. 


Figure  1.  Fallout  Contour  From  a  Constant  Wind  (Ref  10:2) 

This  smearing  of  the  fallout  cloud  on  the  ground  surface  produces  an 
activity  distribution  per  unit  area,  A(x,y,t)  ,  that  can  be  found  from 

t 

A(x,y,t)=A(t)J’  f(x,y,t)g(t)dt  (1) 

^  © 

where  A^(t)  represents  the  total  activity  in  the  cloud  at  time  t  and 
g(t)  is  the  fractional  arrival  rate  of  activity  on  the  ground  (Ref  3:207). 

The  function  f(x,y,t)  represents  the  horizontal  distribution  of 
activity  in  the  cloud.  It  is  a  dual  normal  function  with  a  standard 
deviation  in  both  the  downwind  or  x  direction  and  the  crosswind  or  y 
direction  (Ref  3:208).  These  standard  deviations  are  functions  of  the 
fallout  cloud  arrival  time  on  the  ground,  t  =  x/v  .  The  normalized 

X 

horizontal  activity  distribution  is  represented  by 
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where  v  is  assumed  to  be  a  constant  wind  velocity  in  the  x  direction. 

X 

Equation  (1)  is  the  fundamental  equation  used  to  compute  fallout 

footprints  from  smearing  codes.  The  function  g(t)  is  the  key  element  of 

this  equation.  An  example  of  a  surface  burst  g(t)  calculation  is  illus¬ 
trated  in  Figure  2  (Ref  3^2  ).  The  activity  on  the  ground  per  unit 

area  can  then  be  converted  to  a  dose  rate  in  rads  per  hour  at  a  detector 

3  feet  off  the  ground 

D(x,y,t)  =  C  A(x,y,t)  (3) 

where  D(x,y,t)  is  the  dose  rate  and  C  is  a  constant  having  units  of 
rads  per  hour  per  unit  activity  per  unit  area  (Ref  3:207).  The  deriva¬ 


tion  of  the  constant  C  can  be  found  in  Appendix  A. 


III.  The  Stabilized  Nuclear  Cloud 


Background 

The  detonation  of  a  nuclear  weapon  produces  a  very  hot  cloud  of 
radioactive  weapon  debris  and  incandescent  air  called  the  fireball. 

This  fireball  rises,  expands,  and  cools  mainly  by  radiation  convection 
mixing  of  hot  and  cool  air.  When  thermodynamic  equilibrium  is  reached 
with  the  ambient  atmosphere,  the  cloud  ceases  to  rise  and  becomes  stabi¬ 
lized. 

In  a  surface  burst,  large  amounts  of  surface  material  are  lofted 
into  the  cloud  and  vaporized.  As  the  fireball  cools,  the  radioactive 
constituents  become  mixed  and  incorporated  into  the  surface  debris  by 
condensation.  The  resulting  particles  then  fall  to  the  earth.  Figure  3 
illustrates  the  DELFIC  number-size  distribution  of  particles  resulting 
from  a  surface  burst.  In  an  airburst,  where  the  fireball  does  not  inter¬ 
sect  the  ground,  there  are  no  large  amounts  of  surface  debris  lofted  into 
the  cloud.  As  a  result,  the  particles  are  much  smaller  and  less  widely 
distributed  in  size  as  shown  in  Figure  4.  These  particles  descend  to  the 
earth  at  a  much  slower  rate  than  particles  produced  by  a  surface  burst 
(Ref  8:36). 

Cloud  Center  Height 

The  WSEG  and  AFIT  smearing  codes  allow  fallout  particles  to  begin 
their  descent  from  an  altitude  that  corresponds  to  the  initial  stabilized 
cloud  center  height.  This  stabilized  cloud  center  height  is  modeled 
using  an  empirical  function  of  the  form 


HC  =  44.0  +  6.1  InY  -  .205(lnY  +  2.42)  [  InY  +  2.42  |  (4) 

where  Y  is  the  weapon  yield  in  megatons  and  HC  is  the  cloud  center 
height  in  Icilofaet  (Ref  3:213).  This  relation,  when  applied  to  yields 
of  1  kiloton  or  less,  produces  questionable  results  as  is  pointed  out  in 
the  WSEG  documentation  (Ref  14:31).  For  example,  a  weapon  having  a  0.2 
kiloton  yield  is  predicted  to  have  a  stabilized  cloud  center  height  of 
-0.33  kilofeet. 

A  remedy  to  this  dilemma  was  to  compile  visible  cloud  data  from  22 
low  yield  airburst  tests  conducted  between  1945  and  1962.  This  data  was 
originally  compiled  by  the  old  Defense  Atomic  Support  Agency  (DASA)  and 
contained  in  its  1251  series  of  reports  (Ref  9).  The  data  is  summarized 
in  Table  I.  A  polynomial  least  squares  fit  of  the  data  was  performed  to 
obtain  a  cloud  center  height  relation  as  a  function  of  yield.  It  was 
assumed,  as  in  the  WSEG  model,  that  the  center  of  the  radioactive  cloud 
is  located  at  the  same  altitude  as  the  bottom  of  the  visible  cloud 
(Ref  14:24).  The  resulting  fit  is  a  polynomial  of  fourth  degree  having 
the  form 

HC  =  8.30476298  +  .63632921  In  Y  -  .39763749  (In  Y)* 

-  .01364081  (In  Y)3  +  .00560846  (In  Y)**  (5) 

where  Y  is  the  weapon  yield  in  kilotons  and  HC  is  the  stabilized  cloud 
center  height  in  kilofeet.  Unless  otherwise  stated,  Y  will  have  the 
dimensions  of  kilotons.  Figure  5  illustrates  this  functional  fit  along 
with  the  data  extracted  from  DASA  1251.  It  should  be  noted  that  equa¬ 
tion  (5)  provides  best  results  for  yields  ranging  between  0.01  kilotons 


As  mentioned  earlier,  the  key  element  in  computing  grounded  fallout 
activity  in  smearing  codes  is  the  function  g(t).  An  illustration  of 
g(t)  for  an  airburst  is  shown  in  Figure  6.  This  g(t)  might  be  compared 
to  Figure  2  for  a  surface  burst  in  Section  II.  It  can  be  seen  that  the 
airburst,  with  its  smaller  particles,  leads  to  a  steady  fall  over  months 
instead  of  hours.  This  would  occur  in  a  quiescent  atmosphere.  In  actu¬ 
ality,  the  atmosphere  is  very  turbulent.  Airburst  particles  may  experi¬ 
ence  downdrafts,  updrafts,  etc.  that  may  preclude  their  being  grounded. 

The  only  mechanism  left  for  grounding  of  the  activity  is  scavenging  of 
the  radioactive  cloud  by  rain  (or  snow).  The  conclusion  to  be  drawn 
here  is  that  viscous  fall  of  very  small  particles  is  a  poor  model  for 
predicting  airburst  fallout.  Instead,  rain  scavenging  might  be  a  more 
appropriate  model.  This  section  describes  the  computer  model  that  pre¬ 
dicts  the  maximum  dose  rate  due  to  rain  scavenging  that  can  be  encountered 
for  an  airburst. 


where  f(x)  and  f(y)  are  normalized  Gaussian  distributions  in  the  x  and 

y  directions  and  A(z,t  )  is  the  vertical  distribution  of  cloud 

s 


in  curies/m  at  stabilization  time  t  .  At  later  times,  the  cloud  is 

s 

displaced  in  the  x  direction  by  a  constant  wind  of  velocity  v  .  It  is 
assumed  that  there  is  no  y  component  of  wind  velocity.  The  spatial 
distribution  of  activity  is  then  given  by 

A(x,y,z,t)  =  A(z.t)  f(x)  f(y)  (7) 

If  it  is  assumed  a  rain  cloud  washes  out  the  radioactive  cloud  at  time  t 
by  instantaneously  depositing  100  percent  of  the  activity  onto  the  ground, 
then  the  horizontal  distribution  of  activity  on  the  ground  per  unit  area 
is  described  by 

o 

A(x,y,t)  =  [J  A(z,t)dz]  f(x)  f(y)  (8) 

z 

The  upper  limit  on  the  integral  is  taken  to  be  zero  since  the  integra¬ 
tion  is  carried  out  from  the  top  of  the  cloud  z  (at  time  t)  to  the  ground. 

Particle  Number-Size  and  Activity-Size  Distributions 

Section  HI  stated  that  airburst  particles  are  much  smaller  than 
surface  burst  particles.  Airburst  particles  are  assumed  to  be  spherical 
in  shape  and  obey  a  lognormal  distribution.  The  number-size  distribution 
as  a  function  of  particle  radius  is  expressed  as 


N(r) 


where  r  is  the  particle  radius  in  microns.  The  shaping  parameters,  ao 


and  6o .  are  given  by 


Table  II 


Air  Burst  Number-Size  Distributions 


DISTRIBUTION  NAME 

GEOMETRIC  STANDARD 

DEVIATION 

GEOMETRIC  MEAN 

PARTICLE  RADIUS 

(microns ) 

Norment 

2.0 

.075 

B 

1.63 

.  105 

C 

1.77 

.06 

D-Sample  1 

1.66 

.  10 

D- Sample  2 

1.67 

.09 

E 

2.25 

.055 

F 

1.85 

.043 

G 

2.16 

.0325 

H 

1.92 

.0385 

Bo  =  ln(  o  ) 


where  r  ^  is  the  geometric  mean  particle  radius  in  microns  and  o  is  the 
geometric  standard  deviation  of  the  distribution  (Ref  2).  Airburst 
number-size  distributions  have  been  found  by  Norment  (Ref  13:46)  and 
Nathans  (Ref  11:7567).  The  parameters  for  these  distributions  are  shown 
in  Table  II.  The  distribution  proposed  by  Norment  with  =  ln(.075) 


and  6o 


ln(2.0)  will  be  used  for  all  calculations  made  in  this  thesis 


The  activity-size  distribution  is  expressed  as  a  weighted  sum  of  two 
lognormal  distributions 


A(r)  =  f  A  +  ( 1  -  f  )A 

V  V  vs 


(12) 


where  A  is  the  volumetric  activity-size  distribution,  A  is  the  surface 
V  s 

activity-size  distribution,  and  f^  is  the  volumetric  fractionation  ratio 
of  total  activity  (Ref  3:210).  The  distributions  A  and  A  are  propor- 

V  s 

tional  to  the  third  and  second  moments  of  the  particle  size  distribution 
respectively.  Since  the  nth  moment  of  a  lognormal  is  also  a  lognormal 
distribution  with  the  same  value  for  3.»  values  of  can  be  com¬ 

puted  from  the  relationship  given  by  Aitcheson  and  Brown  (Ref  2:12). 


a  =  Oo  +  n3^ 


(13) 


Russell  (Ref  16:18)  has  suggested  that  smaller  particles,  typical  of 
those  found  in  airbursts,  are  often  volumetrically  distributed  in  activ¬ 
ity.  Talcing  this  suggestion  and  assuming  that  activity  in  airburst 
particles  is  totally  distributed  volumetrically,  f^  =  1.0  ,  equation  (12) 

simplifies  to 


(14) 


where  Oj  =  ln(.317)  and  =  ln(2.0)  .  This  is  the  activity-size 

distribution  that  will  be  used  in  this  thesis.  The  cumulative  activity- 


size  distribution  is  shown  in  Figure  7. 


Calculation  of  A(2,t) 


The  vertical  distribution  of  cloud  activity  per  meter  of  altitude 
is  computed  using  a  method  developed  by  Bridgman  and  Hickman  (Ref  4:2) 
and  is  expressed  as 


A(z.t) 


00 

A( z, r, t )dr 


(15) 


where  A(z,r,t)  is  the  specific  activity  in  curies  per  vertical  meter  of 
altitude  per  micron  of  radial  size  at  time  t.  The  integral  in  equation 
(15)  can  be  replaced  by  a  summation  over  a  set  number  of  discrete  activity 
size  groups,  each  group  containing  an  equal  percentage  of  the  total  cloud 
activity  at  time  t,  where  it  is  assumed  each  group  contains  monosized 
particles  of  mean  radius  r^. 

Since  the  activity-size  distribution  shown  in  Figure  7  has  a  small 
range  of  particle  sizes,  it  can  be  divided  into  10  discrete  size  groups, 
A^(t),  each  having  1/10  of  the  total  cloud  activity.  Table  III  shows  the 
values  of  the  mean  particle  radii  for  each  of  these  activity-size  groups. 
Equation  (15)  can  now  be  expressed  as 


A(z,t) 


f . (z.t) 


i  =  l 


(16) 


where  A^  is  the  total  cloud  activity  and  f^(z,t)  is  the  normalized  ver¬ 
tical  activity  distribution  of  each  size  group  per  meter  of  altitude  at 
any  given  time  t.  The  total  cloud  activity  can  be  expressed  using  the 


distribution  of  cloud  activity  in  curies /m  is  now  expressed  as 

10 

A(z,t)  =  530  X  10*  Y  ^  fAz.t)  (18) 

The  values  of  f^(2,t)  are  computed  by  assuming  each  size  group  to  be 
initially  situated  at  the  stabilized  cloud  center  height,  HC.  Each  size 
group  distribution  is  then  allowed  to  fall  through  the  atmosphere  as  a 
fixed  body  until  it  impacts  on  the  ground,  compressing  much  like  an 
accordian . 

Fall  Mechanics 

The  fall  velocity  of  each  activity-size  group  can  be  computed 
directly  using  Stoke's  Law  fall  mechanics.  Stoke's  Law  can  be  used 
since  the  mean  particle  radius  of  each  group  is  considerably  less  than 
10  microns.  Stoke's  Law  for  free  falling  spheres  in  a  viscous  medium, 
namely  the  atmosphere,  is 


6iT)jvr=  -trrspg 


(19) 


where  y  is  the  dynamic  viscosity  of  the  atmosphere  in  kg/m-sec,  v  is 
the  particle  fall  velocity  in  m/sec,  r  is  the  particle  radius  in  m,  p 
is  the  particle  density  in  kg/m*,  and  g  is  the  gravitational  accelera¬ 
tion  constant  which  has  been  set  to  9.8  m/sec*.  The  value  for  y  is 
altitude  dependent.  A  method  for  calculating  the  dynamic  viscosity 
utilizes  the  US  Standard  Atmosphere  (Ref  17)  and  is  shown  in  Appendix  B. 
Solving  for  velocity  gives 


21 


By  replacing  v  in  equation  (20)  with  Az/At,  the  altitude  of  each  size 
group  can  be  computed  at  different  times  t  from 

n 

z.  (t)  =  HC  -  ^2  At^  (21) 

^  3*1  ^ 


where 


t  = 


E 


3  =  1 


A  t- 


(22) 


The  vertical  distribution  of  cloud 
culated  from 

10 

A(z,t)  =  530  X  10®  Yt"^'^  ^2 

i  =  l 


activity  in  curies/m  can  now  be  cal- 


(23) 


where  O  =  0.18HC  ,  the  standard  deviation  about  the  size  group  center 

z 

height  altitude  as  specified  in  WSEG  (Ref  14:51). 

Dose  Rate  Calculation 

Equation  (23)  was  evaluated  at  several  different  times.  The  pro¬ 
files  are  shown  in  Figures  8  and  9.  It  can  be  seen  that  for  times  up  to 
5  days  after  cloud  stabilization,  the  total  cloud  activity  three  standard 


deviations,  o  ,  below  the  distribution  center  height  still  has  not 


fallen  below  500  m  altitude.  This  indicates  that  an  immersion  dose  due 
to  an  airburst  is  not  a  significant  threat  on  the  ground.  Since  all  the 
cloud  activity  still  resides  in  the  air,  equation  (8)  can  be  reduced  to 


A(x,y,t)  =  530  X  10®  f(x)  f(y) 


Thus,  equation  (24)  gives  the  ground  distribution  of  activity  per  unit 
area  due  to  total  cloud  washout.  A  worst-case  scenerio  will  allow  cal¬ 
culation  of  the  maximum  ground  distribution  of  activity.  By  letting 
X  =  v^t  and  choosing  y  =  0  ,  y  lying  on  the  fallout  footprint  hot¬ 

line,  the  maximum  activity  grounded  can  be  computed.  The  distributions 
in  the  x  direction  and  y  direction  are  now  simplified  to 


f(x,t)  =  - - - 

/2  TT  a 


f(0,t)  = 


These  distributions  are  functions  of  time  because  the  standard  deviations 

0  and  a  are  functions  of  time.  Values  for  O  and  a  are  computed 
X  y  X  y 

in  meters  using  the  WSEG  formulae  (Ref  14:50). 


a  (t) 

X 


a  (t) 

y 


O  S  t 

z  y  . 


»  /•  «,*• 


where  S  and  S  represent  wind  shear  in  the  x  and  y  directions  respec- 
X  y 

tively.  Wind  shear  is  expressed  in  units  of  km/hr-km.  By  using  the 
assumption  of  a  constant  wind  in  only  the  x  direction  and  a  wind  shear 
only  in  the  y  direction,  equation  (27)  simplifies  to 


a^it)  = 


The  value  for  t*  is  limited  to  a  maximum  value  of  three  hours.  This 

limitation  is  made  because  the  second  term  in  equations  (27)  and  (28) 

represents  the  toroidal  growth  of  the  cloud.  This  toroidal  growth  is 

assumed  to  cease  after  three  hours  (Ref  4:7).  Values  for  n  are  com- 

o 

puted  in  meters  using  the  yield  dependent  relationship  expressed  in 
WSEG  (Ref  14:51) 


0  -  1609  exp  I  0.70  +  -  (In  Y)  -  — — 

o  I  3  4.0 


3.25 _  I 

+  (In  Y  +  5.4)2  I 


where  Y  is  the  weapon  yield  in  megatons.  The  WSEG  time  constant,  T^, 
is  computed  in  hours  using  WSEG's  own  low  yield  correction  (Ref  15:1) 


[“(S)-“(S)'] 


1  -  . 5  exp 


[  (-)ll. 


where  HC  is  the  stabilized  cloud  center  height  in  kilofeet. 

The  maximum  dose  rate  can  now  be  computed  using  the  peak  grounded 
activity.  This  peak  is  illustrated  in  Figure  10.  The  dose  rate  is 
computed  by  multiplying  this  peak  grounded  activity  by  a  constant  as 


Figure  10.  Peak  Grounded  Activity 

shown  in  Appendix  A.  The  dose  rate  is  then  given  by 

D(x,o,t)  =  14.4A(x,0,t)  (32) 

The  dose  rate  can  also  be  expressed  in  terms  of  a  unit  time  reference 
dose,  UTRD,  rate.  This  is  the  dose  rate  one  hour  after  weapon  detona¬ 
tion.  The  UTRD  rate  is  expressed  as 

UTRD  =  D(x,0,t)t^  (33) 


V.  Results 


Base  Case 

A  base  case  was  established  as  a  comparison  standard.  Yield  was 
set  at  1.0  kiloton  using  the  Norment  particle  number-size  distribution, 
ao  =  I*'  (-075)  and  =  In  (2.0)  .  Fallout  particle  density  was 

chosen  to  be  5500  kg/m^  (Ref  11:7565).  All  base  case  parameters  are 
summarized  in  Table  IV.  The  cloud  washout  time  after  stabilization  was 
allowed  to  vary  from  1  to  36  hours.  The  x  coordinate  of  the  cloud  center 
point  was  set  to  a  value  v  t  and  then  varied  by  factors  of  a  ,2a, 

X  XX 

and  3  a^.  Since  the  f(x)  distribution  is  a  normal  Gaussian,  the  dose 
rate  and  dose  to  infinity  values  are  symmetric  about  the  cloud  center 

point  in  the  x  direction.  The  results  of  dose  rate  calculations  for  ten 
different  yields  are  summarized  and  graphically  displayed  in  Figures  11 
through  20.  Infinite  dose  calculations  for  these  same  yields  are  shown 
Figures  22  through  31  and  can  be  found  in  Appendix  D. 


Table  IV 


Base  Case  Parameters 


YIELD  =  1.0  KILOTON 


WIND  SHEAR  IN  THE  Y  DIRECTION  =  1.0  KM  PER  KM  PER  HOUR 


DOWNWIND  VELOCITY  =  4.0  METERS  PER  SECOND 


VOLUMETRIC  FRACTIONATION  RATIO  =  1.00 


PARTICLE  MASS  DENSITY  =  5500.  KILOGRAMS  PER  CUBIC  METER 
ALPHAO  =  LN(.075)  =  -2.590  BETAO  -  LN(2.0)  =  .693 

ALPHA2  =  LN(.196)  =  -1.629  ALPHA3  =  LN(.317)  =  -1.149 


Parameter  Variations 


Several  airburst  particle  number-size  distributions  shown  in  Table  II 
were  input  into  the  rain  scavenging  model.  The  resulting  dose  rates  and 
infinite  doses  reflect  no  change  from  the  base  case.  This  result  is  not 
surprising  because  of  the  very  small  radii  of  the  particles.  The  parti¬ 
cles  will  remain  aloft  for  many  days  until  scavenged  from  the  cloud,  re¬ 
gardless  of  the  number-size  distribution. 

A  change  in  particle  mass  density  did  not  affect  the  dose  rates  or 
infinite  doses.  Again,  this  result  is  not  surprising  for  the  same  reason 
stated  above. 

Downwind  velocity  was  changed  over  a  wind  range  of  values.  As 
expected,  dose  rates  and  infinite  doses  did  not  change.  Changes  in  down¬ 
wind  velocity  only  displace  the  cloud  center  point  to  different  locations. 

Different  values  for  wind  shear  in  the  y  direction  were  inserted  into 
the  model.  The  dose  rates  and  infinite  doses  calculated  did  change.  As 
expected,  a  decrease  in  the  wind  shear  results  in  an  increase  in  the  dose 
rate.  As  wind  shear  decreases,  less  cloud  activity  is  displaced  in  the 
crosswind  direction  from  the  cloud  center  point  on  the  ground  allowing 
more  activity  to  be  deposited  during  scavenging. 

Summary 

The  rain  scavenging  model  provides  a  simple  tool  useful  in  predict¬ 
ing  low  yield,  airburst  fallout.  A  tactical  ground  commander  need  only 
know  local  wind  data,  weapon  yield,  ground  zero,  and  time  of  actual  or 
anticipated  rain.  With  this  information,  maximum  dose  rates  and  doses 
at  locations  upwind  and  downwind  from  the  center  of  grounded  activity 
can  be  determined. 
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VI.  Conclusions  and  Recommendations 


Conclusions 

Based  on  the  results  of  this  study,  the  following  conclusions  can 
be  drawn. 

a.  Utility  of  the  model  is  limited  to  weapon  yields  ranging  between 
0.01  kilotons  to  1.0  kilotons.  This  results  from  the  limited  low  yield 
visible  cloud  data  available  from  atmospheric  tests. 

b.  Fractional  arrival  rate  of  activity  function,  g(t),  presently 
used  in  smearing  codes,  is  a  poor  model  for  predicting  airburst  fallout. 

c.  The  immersion  dose  resulting  from  an  airburst  is  not  significant 
as  a  ground  threat. 

d.  Airburst  fallout  is  only  tactically  significant  through  the 
mechanism  of  cloud  scavenging  up  to  36  hours  after  cloud  stabilization. 

Recommendations 

The  rain  scavenging  model  provides  a  useful  airburst  fallout  pre¬ 
diction  tool;  however,  there  is  room  for  improvement.  The  following 
recommendations  would  be  appropriate  for  future  study  and  are  hereby 
offered. 

a.  Compile  a  larger  cloud  data  base  with  the  intent  of  developing 
a  single  yield  dependent  cloud  center  height  relation  useful  from  sub- 
kiloton  to  megaton  yields. 

b.  Incorporate  a  raindrop  size  distribution  into  the  model  and 
begin  the  fallout  deposition  process  at  the  time  rainfall  begins.  The 
intent  is  to  compute  a  raindrop  activity  deposition  function,  g(t). 


A(x,y,t)  can  be  ignored  and  the  integration  taken  to  infinity.  Addition¬ 


ally,  the  gamma  ray  photons  under  go  absorption  by  air  molecules,  are 
spherically  attenuated  as  distance  away  from  the  detector  increases,  and 
are  attenuated  by  one  other  absorbing  media,  namely  the  detector.  The 
exposure  rate  at  ground  position  x,  y  and  at  time  t  can  now  be  represented 
by 


D(x,y,t)  = 


/  (’t)  (hi) 


dA 


(A-1) 


where 


D(x,y,t)  =  Exposure  rate  at  time  t  (roentgens/hr) 
A(x,y,t)  =  Ground  distribution  of  activity  (curies/m^) 


=  Mass  absorption  coefficient  for  air  (m^/kg) 

=  Total  absorption  coefficient  of  the  detector  (m  ^ ) 
=  Distance  from  point  emitter  to  detector  (m) 

=  3.7  X  10 ^ ° disintegrations/sec-curie 
=  1.6  X  10  ^ ^ joules/disintegration 

=  0.00877  joules/kg-roentgen 


By  substituting  dA  =  2  it  rdr  ,  assuming  the  average  gamma  ray  photon 


energy  to  be  1  MeV,  and  using 
equation  (A-l)  reduces  to 


D(x,y,t)  =  3.4A(x,y,t 


=  0.0028  m2/kg 


Mt« 


rdr 


(Ref  7:713), 


(A-2) 


where  the  dimensions  on  the  constant  are  in  roentgens  per  hour  per  unit 


activity  per  unit  area. 

The  above  integral  can  be  simplified  further.  The  geometry  in 
Figure  21  yields: 


s2  =  (lm)2  +  r^  (A-3) 

Differentiating  both  sides,  substituting  into  the  integral,  and  redefin¬ 
ing  integration  limits  gives: 


o 


(A-4) 


Now,  let: 


(A-5) 


and 

ds  dz 

s  "  z  (A-6) 


Substituting  into  equation  (A-4)  and  simplifying  gives: 


e 


rdr 


(A-7) 


The  integral  on  the  right  hand  side  of  equation  (A-7)  is  now  in  the  form 
of  an  Exponential  Integral  of  the  First  Kind,  Ej  (  y^).  For  air,  the 
value  of  is  taken  to  be  0.0082  m  (Ref  7:689).  Substituting  and 

using  an  approximation  for  Ej  (.0082)  from  Abramowitz  and  Stegun  (Ref  1:231) 


gives 


/ 


-2 


dz  »  4.23 


(A-8) 


.0  0  82 


Substituting  this  result  into  equation  (A-2)  yields: 


D(x,y,t)  =  14.4A(x,y,t) 


(A-9) 


This  gives  the  exposure  rate  in  roentgens  per  hour. 

When  computing  the  dose  rate  in  rads  per  hour,  the  detector  would 
be  human  tissue  rather  than  air.  Since  human  tissue  is  comprised  mainly 
of  water,  the  value  of  water  can  be  used  as  an  approximation.  The 

value  of  the  constant  also  changes  to  0.01  joules/rad.  The  end  result 
is  approximately  the  same  as  equation  (A-9).  Equation  (A-9)  can  then  be 
used  to  convert  a  ground  activity  distribution  per  unit  area  to  a  dose 


rate  in  rads  per  hour. 


Appendix  B 


Calculation  of  Atmospheric  Dynamic  Viscosity 

The  US  Standard  Atmosphere  (Ref  17)  provides  empirical  relations  of 
atmospheric  properties  as  functions  of  altitude.  These  relations  have 
been  fit  to  atmospheric  data  at  various  levels  in  the  atmosphere. 

The  largest  weapon  yield  utilized  in  this  thesis  is  1  )iiloton.  The 
cloud  center  height  predicted  by  equation  (5)  for  1  kiloton  is  2531  m. 

For  calculating  the  dynamic  viscosity  of  the  atmosphere,  it  is  only 
necessary  to  use  the  relations  given  for  altitudes  up  to  11000  m. 

The  US  Standard  Atmosphere  defines  the  dynamic  viscosity  as 

1.458  X  10~^  T^-^ 

^  °  T  +  110.4  (B-1) 

where  T  is  the  ambient  temperature  of  the  atmosphere  in  'K  and  y  is  in 
kg/m-sec  (Ref  17:19).  The  ambient  atmospheric  temperature  in  °K  is 
expressed  as  a  function  of  altitude  by 

T(z)  =  288.15  -  .0065Z  (B-2) 

where  z  is  the  altitude  in  m  (Ref  17:10).  By  substituting  equation  (B-2) 
into  equation  (B-1),  the  dynamic  viscosity  of  the  atmosphere  is  computed 


directly  from  any  given  altitude  up  to  11000  m. 


Appendix  C 


ALPHAO 

ALPHA2 

ALPHA3 

AXYT 

BETAO 

CODE 

DDOT 

DELTAT 

DELTAZ 

EXPO 

EXPOFX 

EXPOFY 

DOSE 

FV 

FX 

FY 

GA 

HC 


I .  Program  Variable  Definitions 


=  Parameter 

of 

the 

logarithm 

of 

the 

=  Parameter 

of 

the 

logarithm 

of 

the 

=  Parameter 

of 

the 

logarithm 

of 

the 

particle  number  size  distribution: 
median  particle  radius. 

surface  activity-size  distribution: 
median  particle  radius. 

volumetric  activity-size  distribution; 
median  particle  radius. 


=  Ground  distribution  of  activity  per  unit  area  (curies/m^). 


=  Logarithmic  slope  of  the  particle  number-size  distribution 


=  An  integer  code  used  for  computing  different  doses: 
Code  0  for  dose  to  infinity,  code  1  for  dose  to  24 
hours,  code  2  for  dose  to  4  hours. 


=  Dose  rate  (Rads /hr). 


=  Time  step  parameter  set  to  1.0  hours.  It  is  used  in 
the  Stoke 's  Law  fall  mechanics. 


=  Distance  fallout  particles  fall  during  time  DELTAT  (m). 

«  Dummy  variable  representing  the  exponential  argument  in 
the  variable  SIGO. 

=  Dummy  variable  representing  the  exponential  argument  in 
the  variable  FX. 

=  Dummy  variable  representing  the  exponential  argument  in 
the  variable  FY. 


=  Gamma  radiation  dose  (Rads). 

=  Volumetric  fractionation  ratio. 

=  Normalized  Gaussian  distribution  in  the  x  direction. 
=  Normalized  Gaussian  distribution  in  the  y  direction. 
=  Acceleration  due  to  gravity  (9.80  m/sec^). 

=  Stabilized  cloud  center  height  (m). 


MU 


Dynamic  viscosity  (kg/m-sec). 


M 


20. 

NUMBER 

= 

Dummy  variable  representing  the  squared  value  of  SIGY. 

21. 

PI 

= 

The  mathematical  constant  it  . 

22. 

RHOP 

= 

Mass  density  of  fallout  particles  (kg/m^). 

V 

23. 

RMEAN 

= 

Vector  specifying  the  values  of  mean  particle  radii 

•J 

in  each  of  the  10  equal  activity-size  groups  (microns). 

24. 

SHEARX 

= 

Wind  shear  in  the  x  direction  (km/hr-km). 

\ 

25. 

SMEARY 

= 

Wind  shear  in  the  y  direction  (km/hr-km). 

26. 

SIGO 

= 

Standard  deviation  of  horizontal  cloud  radius  (m). 

27. 

SIGX 

= 

Cloud  standard  deviation  in  the  x  direction  (m). 

F 

28. 

SIGY 

= 

Cloud  standard  deviation  in  the  y  direction  (m). 

29. 

SIGZ 

= 

Cloud  standard  deviation  in  the  vertical  direction  (m). 

30. 

TC 

= 

WSEG  time  constant  (hours). 

r. 

31. 

TEMP 

s 

Ambient  atmospheric  temperature  (“K). 

^ ' 

32, 

TIME 

= 

Time  of  cloud  washout  (hours). 

• 

* 

33. 

TSTAR 

WSEG  toroidal  growth  time  of  cloud  (hours). 

34. 

UTRD 

= 

Unit  time  reference  dose  rate  (Rads/hr). 

,  * 

35. 

V 

= 

Fallout  particle  fall  velocity  (m/sec). 

36, 

VX 

= 

wind  velocity  in  the  x-direction  (m/sec). 

*« 

rn 

1 

37. 

VY 

= 

Wind  velocity  in  the  y-direction  (m/sec). 

38. 

XO 

= 

Cloud  center  position  in  the  x-direction  (m). 

•J 

39. 

X 

= 

Distance  east/west  of  XO  (m). 

■ 

3 

40. 

YO 

= 

Cloud  center  position  in  the  y-direction  (m). 

41. 

Y 

= 

Distance  north/south  of  YO  (m). 

J* 

• 

42. 

YIELD 

= 

Weapon  yield  (kilotons). 

t' 

43. 

ZC 

= 

Two-dimensional  array  of  particle  size  group  center 

4 

1 

altitudes  at  different  times  of  fall. 

> 

•- 

> 

I 


RAINOUT  MODEL 


This  program  computes  the  maximum  dose  rate.  In  Rads  per  hour, 
produced  by  a  nuclear  air  burst.  The  nuclear  cloud  is  washed  out  by  a 
simulated  r-fn  storm  at  a  specified  time  t  of  cloud  washout.  The 
activity  grounded  as  a  result  of  the  radioactive  cloud  washout, 
A(x,y,t>,  Is  used  to  compute  the  dose  rate.  This  program  will  accomo¬ 
date  changes  In  the  following  Input  parameters:  particle  size  dist¬ 
ribution,  weapon  yield  (for  yields  ranging  from  .01  to  1.0  kllotons), 
fractionation  ratio,  fallout  particle  density  (in  kilograms  per  cubic 
meter), wind  velocity  (In  meters  per  second),  and  wind  shear  (In 
kilometers  per  hour  per  kilometer  of  atmospheric  altitude). 


REAL  PI .ALPHA0,EETA0.FV,RHOP .GA.RMEAN( 101 . Y I  EL D , HC . AL P HA2 , AL P HA3 
REAL  TEMP .MU . DELTAZ . VPR I  NT . X0, Y0 . AXVT , DDOT ( 4 ) , DOSE( 4 ) 

REAL  V.SIGMAZ( 1 0 ), DELTA! . ZC ( 10. 0: 50 > , T I  ME ( 0 : 50 ) , UTR D ( 4 ) 

REAL  TC .SIGZ ,SIG0. EX PO.SIGX .TSTAR ,SIGY. NUMBER . SHEARX . SMEARY , VX , VY 
REAL  X( 4 ) . Y. FX .FY, EXPOFX.EXPOFY 
INTEGER  1  . 0 .K.CODE .K1 

INPUT  OF  ALL  PROBLEM  PARAMETERS. 

PARAMETER  (PI=3. 14159,  AL PHA0=LOG ( . 075 ) .  BETA0=LOG ( 2 . 0 ) ,  FV=1.0. 

+  RHOP=5500.0,  GA=9.80.  YIELD=0.1.  VX=4.0,  VY=0.0, 

+  SHEARX=10.0,  SHEARY=1.0.  Y=0.0.  CODE=0  ) 

OPEN  ( 3.FILE=’dosedata’ ) 

REWIND  3 


iLd 


SPECIFY  THE  MEAN  PARTICLE  RADII  FOR  THE  TEN  EQUAL  ACTIVITY  GROUPS. 

CALL  MEAN'ALPHA0.BETA0.ALPHA2.ALPHA3.RMEAN  > 

COMPUTE  THE  ALTITUDE  OF  THE  STABILIZED  CLOUD  CENTER.  IN  METERS. 

HC=( 1 609. /5. 28 )»( 8. 30476298+. 63632921 ‘LOGI YIELD )-.39763749»LOG(YIE 
+LD ) ‘LOG! YIELD )- .0136408 1 ’LOG (YIELD )**3+ .00560846»LOG(Y1 ELD )**4 ) 
TC=( 12.*(HC/60. )*{5.28/1609. ) -2 . 5* ( HC / 60 . ) • ( 5 . 28 / 1 609 . >*(5.28/1609 
+  .))•(!. 0-.5»EXP((HC/25. >*(5.28/1609. ) * ( 5 . 28/ 1 609 . > ) ) 

SIGZ=. 18*HC 

EXPO=.7+(l./3.  >*LOG( YIELD/ 1000. ) - ( 3 . 25 / 1  ■; .  +  ( LOG  YIELD/ 1000. >+5 . 4  > 

+  *( LOG( YIELD/ 1000. )+5.4 )) > 

SIG0=I609.*EXP(EXPO) 

PRINT  1.  YIELD 

1  FORMAT  (/."YIELD  =  "  ,  F 6 . 4  ,  "  KILOTON  AIRBURST") 

PRINT  2,  HC 

2  FORMAT  (/."STABILIZED  CLOUD  CENTER  HEIGHT  =  ".F5.0,"  METERS") 

PRINT  25.  TC 

25  FORMAT  (/."TC  *  ",F5.2,"  HOURS") 

PRINT  26.  SIG0 

26  FORMAT  (/,"SIGMA0  =  ".F5.0,"  METERS") 

PRINT  27.  SIGZ 

27  FORMAT  (/."SIGMAZ  =  ".F4.0."  METERS") 

PRINT  28.  VX 

28  FORMAT  (.•."DOWNWIND  VELOCITY  =  ".F4.]."  METERS  PER  SECOND") 

PRINT  3.  FV 


■W  .S'wV.V  ••  -N  .V  -'kW.V V' V  -v’ 
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F4 . 2  ) 


3  FORMAT  (/."VOLUMETRIC  FRACTIONATION  RATIO  =  ", 

PRINT  35  RHOP 

35  FORMAT  (/."PARTICLE  MASS  DENSITY  =  “,F5.0.”  KILOGRAMS  PER  CUBIC  ME 
♦  TER"  > 

PRINT  4.  ALPHAJ0f.BETA0 

4  FORMAT  (/,"ALPHA0  =  LN(.075)  =  ”,F6.3.“  BETA0  =  LN(2.0>  =  ". 

+F5.3) 

PRINT  5.  ALPHA2. ALPHAS 

5  FORMAT  (/."ALPHAS  =  LN(.196)  =  "  F6.3."  ALPHAS  =  LN(.317>  =  " 
+.F6.3.//) 

COMPUTE  THE  ALTITUDE  FOR  THE  CENTER  OF  EACH  PARTICLE  SIZE  GROUP. 

DO  7  1=1.10 
ZC( I .0)=HC 

SIGMAZ( I )=. 18*ZC( I .0) 

7  CONTINUE 

COMPUTE  THE  FALL  TIMES  FOR  EACH  SIZE  GROUP  CENTER  USING  STOKE’S  LAW 

FALL  MECHANICS. 

TIME(0)=0.0 

DELTAT=1 .0 

DO  13  0=1.36 

TIME(0 )=TIME(0-1 )+DELTAT 
DO  12  1=1.10 

OELTAT=DELTAT*3G00.0 
TEMP=288. 15- .0065*ZC( I .0-1 ) 

MU=( 1 .458E-6*TEMP**1 . 5 > / ( TEMP+ 1 10 . 4 ) 

V=(2.0/9.0;*( (RMEAN( I > *  1 E-6 > **2*8  HOP *GA ) /MU 
DELTAZ=V*DELTAT 
ZC( 1 .0)=2C( I .0-1 >-DELTAZ 
IF (ZC( I .0 ) .LT.0.0)THEN 
2C( I ,0 )=0.0 
GO  TO  105 
END  IF 

05  OELTAT=OELTAT/3600.0 

VPRINT=V*100.0 

12  CONTINUE 

13  CONTINUE 

COMPUTE  THE  VALUE  OF  THE  DOSE  RATE  IN  R/hr. 


1F(CODE.EQ.0)THEN 

PRINT*. •  ’ 

PRINT*. 'DOWNWIND  SIGMAX 

TIME  OF 

DOSE 

+  TO’ 

PRINT*. 'DISTANCE  DISTANCE 

WASH  OUT 

DOSE  RATE 

INFI 

+NITY’ 

PRINT*.'  (Km)  (Km) 

( hours ) 

(Rads/hr ) 

(Ra 

+  ds  )  ’ 

PRINT*.’  ’ 

GO  TO  200 

ELSE  IF(CODE.EQ. 1  )THEN 

PRINT*.’  ' 

PRINT*. 'DOWNWIND  SIGMAX 

TIME  OF 

D 

♦OSE  TO’ 

PRINT*. 'DISTANCE  DISTANCE 

WASH  OUT 

DOSE  RATE 

2 

+4  HOURS’ 

PRINT*.’  (Km)  (Km) 

( hours ) 

( Rads/hr ) 

+ ( Rads ) ’ 

PRINT*.’  ’ 

GO  TO  200 

PRINT*,’  ’ 

PRINT*. ’DOWNWIND  SIGMAX  TIME  OF 

♦DOSE  TO’ 

PRINT*, ’DISTANCE  DISTANCE  WASH  OUT  DOSE  RATE 

+4  HOURS’ 

PRINT*,’  (Km)  (Km)  (hours)  (Rads/hr) 

+  (Rads)  ’ 

PRINT*.’  ’ 

GO  TO  200 

END  IF 

200  DO  17  0=1 .36 

IF (TIME( 0 ) . LE . 3 .0)THEN 
TSTAR=TIME( 0 ) 

GO  TO  135 
END  IF 
TSTAR=3.0 

135  SIGX=SIG0*SQRT( 1 . 0+8 . 0*TSTAR /TC ) 

NUMBER,=SIG0*SIG0*(  1  . 0+3 . 0*TSTAP. /TC  )  +  (  S  I  GZ*SHEARY*T  I  ME  (  0  )  >* 

+  {SIGZ*SHEARY*TIME(0 ) ) 

SIGY=SQRT(NUM8ER ) 

X0=VX*TIME(O )* 3600.0 
Y0=VY*TIME(O )*3G00.0 
DO  14  K=1 ,4 
X(K)=X0+(K-1 )*SIGX 

EXPOFX=-.5*( (X(K)-X0)/SIGX)*! (X(K)-X0)/SIGX> 

EXPOFY=-.5*( (Y-Y0)/SIGY)*( (V-Y0>/S1GY) 

FX=( 1 .0/(SQRT{2.0*PI )*SIGX) )*EXP(EXPOFX) 

FY=( 1 .0/( SQRT( 2 .0*P I)*S1GY))*EXP(EXP0FY) 
AXYT=(530.E6*YIELD*TIME(O )**( -1 .2) )*FX*FY 
DDOT(K)=14.4*AXYT 
UTRD(K)=DDOT{K)»TIME(0)**1 .2 
IF(CODE.EQ.0)THEN 

DOSE (K)=5.0*UTRD(K)*T1ME(O )**(-.  2) 

GO  TO  14 

E(.SE  IFICOOE.EQ.  1  )THEN 

DOSE{K)=5.0*UTRD(K)*(TIME(O )**(-. 2 )-(TIME(0 )+24.0)** 
+  (-.2) ) 

GO  TO  14 

ELSE 

DOSE(K>=5.0*UTRD{K>*(TIME(O )**(-. 2 )-(TIME(0  >  +  4.0>** 

+  (-.2)  > 

GO  TO  14 

END  IF 

14  CONTINUE 
SIGX=SIGX/1000.0 

WRITE  (3,136, I0STAT  =  K1 .ERR=17)  T1ME(0  )  .DDOT( I ) .DDOT( 2) ,DDOT( 3) 
+  DDOT( 4 ) ,DOSE( 1 ) .DOSE( 2) .DOSE( 3 ) ,DOSE( 4 ) ,SIGX 

136  FORMAT  (F4.1. IX. Ell. 6. IX. Ell. 6. IX. Ell. 6. IX. Ell. 6. IX. Ell. 6, IX. 

+  Ell. 6, IX, Ell. 6. IX. Ell. 6. IX. F8. 2) 

PRINT  15,  X.SIGX.TIME(0).OOOT.DOSE 

15  FORMAT  (F8.3.4X,F8.3.9X,F4. 1 .7X.F7.2,5X.F8.2) 

17  CONTINUE 

ENDFILE  (3) 

CLOSE  (3) 


SUBROUTINE  FOR  MEAN  PARTICLE  RADII 


SUBROUriNE  MEANCALPHAjB.BETAB.ALPHAZ.ALPHAS.RMEAN) 
REAL  ALPHAB,BETAB,ALPHA2. ALPHAS. Z< 10 > . RMEAN ( 10 > 
INTEGER  I 

DATA  (Z( I ) . 1=1 , 10>/-1 .64,-1 . 03 . - . 672 . - . 3B4 , - . 126. 
♦1 .03. 1.64/ 

ALPHA2=ALPHA0+2.0*BETA0*BETA0 
ALPHA3=ALPHA0+3.0*BETA0*BETA0 
DO  10  1=1,10 

RMEANl I )=EXP<Z< I ) •BETA0+ALPHA3 > 

10  CONTINUE 
RETURN 
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